####################################################################################
#
#Young Invincibles to Young Insurables: Insurance uptake by the healthy, wealthy and employed
#
####################################################################################

library(AER)
library(foreign)
library(plyr)
library(sandwich)
library(car)
library(zoo)
library(xtable)
library(ggplot2)
library(scales)
library(gridExtra)



rm(list=ls())
#Load Data

if (getwd()=="C:/Users/Trivial/Documents") {setwd ("C:/Users/Trivial/Dropbox/Replication paper/AEJEPcode/")}
if (getwd()=="C:/Users/Akshar/Documents") {setwd ("C:/Users/Akshar/Dropbox/Replication paper/AEJEPcode/")}
if (getwd()=="/Users/CarlosRH") {setwd ("/Users/CarlosRH/Dropbox/Replication paper/AEJEPcode/")}
if (getwd()=="C:/Users/Admin/Dropbox/Harvard course materials/Gov2001/Replication paper/AEJEPcode") {setwd ("C:/Users/Admin/Dropbox/Harvard course materials/Gov2001/Replication paper/AEJEPcode/")}

#Dataset
data<-read.dta("data/depcov_dataset1.dta")

####################################################################################
# DATA PREPARATION
####################################################################################

##Generate Age Dummies

for(i in 16:25){
  name<-paste("age",i,sep="")
  data[,name]<-as.integer(data$age==i)}

for(i in 27:29){
  name<-paste("age",i,sep="")
  data[,name]<-as.integer(data$age==i)}

##Year Dummies
for(i in 2009:2011){
  name<-paste("y",i,sep="")
  data[,name]<-as.integer(data$year==i)}

##Month Dummies
data$month<-as.integer(data$month)
for(i in 1:12){
  name<-paste("month",i,sep="")
  data[,name]<-as.integer(data$month==i)}

rm(name, i)

##Trend Variable
data$trend<-(data$year-2008)*12+(data$month-8)
data$trend2 <- data$trend^2

##Post reform dummies
data$mar_sept10<-as.integer(data$year==2010 & data$month>=3 & data$month<=9)
data$after_oct10<-as.integer((data$year==2010 & data$month>=10) | data$year>=2011)
data$oct10_feb11<-as.integer((data$year==2010 & data$month>=10) | (data$year==2011 & data$month<=2))
data$mar11_sept11<-as.integer(data$year==2011 & data$month>=3 & data$month<=9)
data$after_sept11<-as.integer((data$year==2011 & data$month>9) | data$year>2011)
data$after_mar10<-as.integer((data$year==2010 & data$month>=3) | data$year>=2011)

#Treatment group dummy
data$fedelig<-as.integer(data$age>=19 & data$age<26)

#Interactions Post Reform and Treatment Group
data$elig_mar10<-data$mar_sept10*data$fedelig
data$elig_oct10<-data$after_oct10*data$fedelig
data$elig_middle<-data$oct10_feb11*data$fedelig
data$elig_marsept<-data$mar11_sept11*data$fedelig
data$elig_sept11<-data$after_sept11*data$fedelig

#Splined variables
data$trendbefore_mar10        <- (1-data$after_mar10)*data$trend
data$trendafter_mar10         <- data$after_mar10*((data$year-2010)*12+(data$month-3))

data$trendbefore_mar10_treat  <- (1-data$after_mar10)*data$trend*data$fedelig
data$trendafter_mar10_treat   <- data$after_mar10*((data$year-2010)*12+(data$month-3))*data$fedelig

data$before_mar10 <- as.integer((data$year==2010 & data$month<3) | data$year<=2009)
data$before_mar10_treat <- data$before_mar10*data$fedelig
data$after_mar10_treat <- data$before_mar10*data$fedelig



#Interactions Unempoyment Rate and Treatment group
data$ue_treat<-data$ue*data$fedelig

#Interactions of post-reform dummies, treatment group dummy and an indicator of parents ESI coverage
data$emppar_mar10<-data$emphi_parent*data$mar_sept10
data$emppar_oct10<-data$emphi_parent*data$after_oct10

#Interaction treatment group dummy and parent's ESI coverage
data$fedelig_emppar<-data$fedelig*data$emphi_parent

#Interactions post reform dummies, treatment group dummy and an indicator of parents ESI coverage
data$ddd_mar10<-data$mar_sept10*data$fedelig*data$emphi_parent
data$ddd_oct10<-data$after_oct10*data$fedelig*data$emphi_parent

#Labor Force Status
data$employed<-as.integer(data$labor==1)
data$unemployed<-as.integer(data$labor==2)
data$hour_vary<-as.integer(data$hours==-1)
data[data$hours>-1,"ln_hours"]<-log(data[data$hours>-1,]$hours+1)

##Group Variables
data$age_group<-NA
data[data$age>=16 & data$age<=18,"age_group"]<-1
data[data$age>=19 & data$age<=25,"age_group"]<-2
data[data$age>=27 & data$age<=29,"age_group"]<-3

###Add State trend
v<-sort(unique(data$fipstate))

for(i in 1:length(v)){
  name<-paste("st_trend_",i,sep="")
  data[,name]<-0
  data[data$fipstate==v[i],name]<-data[data$fipstate==v[i],]$trend
}
rm(i, name, v)
### Add State2 trend
v<-sort(unique(data$fipstate))
for(i in 1:length(v)){
  name<-paste("st_trend2_",i,sep="")
  data[,name]<-0
  data[data$fipstate==v[i],name]<-data[data$fipstate==v[i],]$trend2
}
rm(i, name, v)

data$trend_treat<-data$trend*data$fedelig

data$unins<-as.integer(data$anyhi==0)

#Family coverage variable, whether a parent had family coverage in Nov 2009 -Feb 2010 given that a parent had ESI in that period and the young adult did not have dependent coverage through parent's ESI;
data[is.na(data$other_dep_pre),]$other_dep_pre<--1
data[is.na(data$emphi_dep_pre),]$emphi_dep_pre<--1
data[is.na(data$emphi_p_pre),]$emphi_p_pre<--1

data$famcov_pre<-NA
data[data$other_dep_pre==1 & data$emphi_dep_pre==0,]$famcov_pre<-1
data[data$other_dep_pre==0 & data$emphi_p_pre==1 & data$emphi_dep_pre==0,]$famcov_pre<-0


#Sum group variable indicates (pre or post(law)*(treatment or control))

data$sum_group<-NA
data[data$after_oct10==0 & data$mar_sept10==0 & data$fedelig==1,]$sum_group<-1
data[data$after_oct10==0 & data$mar_sept10==0 & data$fedelig==0,]$sum_group<-2
data[data$after_oct10==1 & data$fedelig==1,]$sum_group<-3
data[data$after_oct10==1 & data$fedelig==0,]$sum_group<-4

data$monthdate <- as.Date(paste(data$year,data$month,"1",sep="-"))


####################################################################################
#
# PREPARATION HEALTH
#
####################################################################################
data$healthnum<-as.numeric(data$health)-1

data$ExH<-as.numeric(data$healthnum==1)
data$VGH<-as.numeric(data$healthnum==2)
data$GH<-as.numeric(data$healthnum==3)
data$FH<-as.numeric(data$healthnum==4)
data$PH<-as.numeric(data$healthnum==5)

data$fedeligEx<-data$fedelig*data$ExH
data$fedeligVG<-data$fedelig*data$VGH
data$fedeligG<-data$fedelig*data$GH
data$fedeligF<-data$fedelig*data$FH
data$fedeligP<-data$fedelig*data$PH

data$after_mar10Ex<-data$after_mar10*data$ExH
data$after_mar10VG<-data$after_mar10*data$VGH
data$after_mar10G<-data$after_mar10*data$GH
data$after_mar10F<-data$after_mar10*data$FH
data$after_mar10P<-data$after_mar10*data$PH

data$after_mar10_treatEx<-data$after_mar10_treat*data$ExH
data$after_mar10_treatVG<-data$after_mar10_treat*data$VGH
data$after_mar10_treatG<-data$after_mar10_treat*data$GH
data$after_mar10_treatF<-data$after_mar10_treat*data$FH
data$after_mar10_treatP<-data$after_mar10_treat*data$PH

data$trendafter_mar10Ex<-data$trendafter_mar10*data$ExH
data$trendafter_mar10VG<-data$trendafter_mar10*data$VG
data$trendafter_mar10G<-data$trendafter_mar10*data$G
data$trendafter_mar10F<-data$trendafter_mar10*data$FH
data$trendafter_mar10P<-data$trendafter_mar10*data$P

data$trendafter_mar10_treatEx<-data$trendafter_mar10_treat*data$ExH
data$trendafter_mar10_treatVG<-data$trendafter_mar10_treat*data$VG
data$trendafter_mar10_treatG<-data$trendafter_mar10_treat*data$G
data$trendafter_mar10_treatF<-data$trendafter_mar10_treat*data$FH
data$trendafter_mar10_treatP<-data$trendafter_mar10_treat*data$P


####################################################################################
#
# PREPARATION WEALTH
#
####################################################################################
data$IncQuintile <- factor(cut(data$famyincy, quantile(data$famyincy, probs = seq(0, 1, 0.20)), include.lowest = TRUE), 
                           labels = c(1:5))

data$IncQuintile<-as.numeric(data$IncQuintile)

data$SES1<-as.numeric(data$IncQuintile==1) #LOW QUINTILE
data$SES2<-as.numeric(data$IncQuintile==2)
data$SES3<-as.numeric(data$IncQuintile==3)
data$SES4<-as.numeric(data$IncQuintile==4)
data$SES5<-as.numeric(data$IncQuintile==5) #HIGH QUINTILE

data$fedeligSES1<-data$fedelig*data$SES1
data$fedeligSES2<-data$fedelig*data$SES2
data$fedeligSES3<-data$fedelig*data$SES3
data$fedeligSES4<-data$fedelig*data$SES4
data$fedeligSES5<-data$fedelig*data$SES5

data$after_mar10SES1<-data$after_mar10*data$SES1
data$after_mar10SES2<-data$after_mar10*data$SES2
data$after_mar10SES3<-data$after_mar10*data$SES3
data$after_mar10SES4<-data$after_mar10*data$SES4
data$after_mar10SES5<-data$after_mar10*data$SES5

data$after_mar10_treatSES1<-data$after_mar10_treat*data$SES1
data$after_mar10_treatSES2<-data$after_mar10_treat*data$SES2
data$after_mar10_treatSES3<-data$after_mar10_treat*data$SES3
data$after_mar10_treatSES4<-data$after_mar10_treat*data$SES4
data$after_mar10_treatSES5<-data$after_mar10_treat*data$SES5

data$trendafter_mar10SES1<-data$trendafter_mar10*data$SES1
data$trendafter_mar10SES2<-data$trendafter_mar10*data$SES2
data$trendafter_mar10SES3<-data$trendafter_mar10*data$SES3
data$trendafter_mar10SES4<-data$trendafter_mar10*data$SES4
data$trendafter_mar10SES5<-data$trendafter_mar10*data$SES5

data$trendafter_mar10_treatSES1<-data$trendafter_mar10_treat*data$SES1
data$trendafter_mar10_treatSES2<-data$trendafter_mar10_treat*data$SES2
data$trendafter_mar10_treatSES3<-data$trendafter_mar10_treat*data$SES3
data$trendafter_mar10_treatSES4<-data$trendafter_mar10_treat*data$SES4
data$trendafter_mar10_treatSES5<-data$trendafter_mar10_treat*data$SES5

####################################################################################
#
# PREPARATION EMPLOYMENT
#
####################################################################################

#Interactions Federal Eligibility
data$fedeligemp<-data$fedelig*data$employed
data$fedeligft<-data$fedelig*data$ft
data$fedeligje<-data$fedelig*data$job_ease
data$fedeligemp2<-data$fedelig*data$emp2diffemp
data$fedelighrsv<-data$fedelig*data$hour_vary

data$pt<-as.integer(data$ft==0)
data$fedeligpt<-data$fedelig*data$pt

FedeligE<-as.matrix(subset(data, select=c("fedeligemp", "fedeligft", "fedeligpt", "fedeligje", "fedeligemp2", "fedelighrsv")))

data$after_mar10emp<-data$after_mar10*data$employed
data$after_mar10ft<-data$after_mar10*data$ft
data$after_mar10pt<-data$after_mar10*data$pt
data$after_mar10je<-data$after_mar10*data$job_ease
data$after_mar10emp2<-data$after_mar10*data$emp2diffemp
data$after_mar10hrsv<-data$after_mar10*data$hour_vary

data$after_mar10_treatemp<-data$after_mar10_treat*data$employed
data$after_mar10_treatft<-data$after_mar10_treat*data$ft
data$after_mar10_treatpt<-data$after_mar10_treat*data$pt
data$after_mar10_treatje<-data$after_mar10_treat*data$job_ease
data$after_mar10_treatemp2<-data$after_mar10_treat*data$emp2diffemp
data$after_mar10_treathrsv<-data$after_mar10_treat*data$hour_vary

EffectsSplineLevelsE<-as.matrix(subset(data, select=c("after_mar10emp", "after_mar10ft","after_mar10pt", "after_mar10je", "after_mar10emp2", "after_mar10hrsv", "after_mar10_treatemp", "after_mar10_treatft","after_mar10_treatpt", "after_mar10_treatje", "after_mar10_treatemp2", "after_mar10_treathrsv")))

data$trendafter_mar10emp<-data$trendafter_mar10*data$employed
data$trendafter_mar10ft<-data$trendafter_mar10*data$ft
data$trendafter_mar10pt<-data$trendafter_mar10*data$pt
data$trendafter_mar10je<-data$trendafter_mar10*data$job_ease
data$trendafter_mar10emp2<-data$trendafter_mar10*data$emp2diffemp
data$trendafter_mar10hrsv<-data$trendafter_mar10*data$hour_vary

data$trendafter_mar10_treatemp<-data$trendafter_mar10_treat*data$employed
data$trendafter_mar10_treatft<-data$trendafter_mar10_treat*data$ft
data$trendafter_mar10_treatpt<-data$trendafter_mar10_treat*data$pt
data$trendafter_mar10_treatje<-data$trendafter_mar10_treat*data$job_ease
data$trendafter_mar10_treatemp2<-data$trendafter_mar10_treat*data$emp2diffemp
data$trendafter_mar10_treathrsv<-data$trendafter_mar10_treat*data$hour_vary

EffectsE<-as.matrix(subset(data, select=c("trendafter_mar10emp", "trendafter_mar10ft","trendafter_mar10pt","trendafter_mar10hrsv", "trendafter_mar10je", "trendafter_mar10emp2", "trendafter_mar10_treatemp","trendafter_mar10_treatft", "trendafter_mar10_treatpt", "trendafter_mar10_treathrsv","trendafter_mar10_treatje", "trendafter_mar10_treatemp2")))

Jobmatrix<- as.matrix(subset(data, select=c("employed", "ft", "pt", "hour_vary", "job_ease", "emp2diffemp" )))

####################################################################################
#
# REGRESSIONS VARIABLES
#
####################################################################################
Outcomes<-subset(data, select=c(anyhi, emphi_dep, indiv, emphi, govhi))
Effects<-as.matrix(subset(data, select=c(mar_sept10, after_oct10, fedelig, elig_mar10, elig_oct10)))
EffectsSplineRates <- as.matrix(subset(data, select=c(trendafter_mar10,trendafter_mar10_treat)))
EffectsSplineLevels <- as.matrix(subset(data, select=c(after_mar10,after_mar10_treat)))
ExpVars<-as.matrix(subset(data, select=c(age17, age18, grep("age2", names(data)),trend, female, ue, hispanic, white, asian, other, student, mar, grep("fpl_ratio", names(data)), ue_treat, y2009, y2010, month2, month3, month4, month5, month6, month7, month8, month9, month10, month11, month12, grep("st_trend_", names(data)))))
ExpVarsNoTime  <- as.matrix(subset(data, select=c(fedelig,age17,age18,age19,age20,age21,age22,age23,age24,age25,age27,age28,age29,female,hispanic,white,asian,other,student,mar,ue,fpl_ratio,fpl_ratio_2,ue_treat)))
StateTrends<-as.matrix(subset(data, select=c(grep("st_trend_", names(data)))))
StateTrend2s<-as.matrix(subset(data, select=c(grep("st_trend2_", names(data)))))


EffectsH<-as.matrix(subset(data, select=c("trendafter_mar10Ex", "trendafter_mar10VG", "trendafter_mar10G", "trendafter_mar10F", "trendafter_mar10P","trendafter_mar10_treatEx","trendafter_mar10_treatVG", "trendafter_mar10_treatG", "trendafter_mar10_treatF", "trendafter_mar10_treatP")))
EffectsSplineLevelsH<-as.matrix(subset(data, select=c("after_mar10Ex", "after_mar10VG", "after_mar10G", "after_mar10F", "after_mar10P", "after_mar10_treatEx", "after_mar10_treatVG", "after_mar10_treatG", "after_mar10_treatF", "after_mar10_treatP")))
Health<-as.matrix(subset(data, select=c("ExH", "VGH", "GH", "FH", "PH")))
FedeligH<-as.matrix(subset(data, select=c("fedeligEx", "fedeligVG", "fedeligG", "fedeligF", "fedeligP")))

EffectsW<-as.matrix(subset(data, select=c("trendafter_mar10SES1", "trendafter_mar10SES2", "trendafter_mar10SES3", "trendafter_mar10SES4", "trendafter_mar10SES5","trendafter_mar10_treatSES1","trendafter_mar10_treatSES2", "trendafter_mar10_treatSES3", "trendafter_mar10_treatSES4", "trendafter_mar10_treatSES5")))
EffectsSplineLevelsW<-as.matrix(subset(data, select=c("after_mar10SES1", "after_mar10SES2", "after_mar10SES3", "after_mar10SES4", "after_mar10SES5", "after_mar10_treatSES1", "after_mar10_treatSES2", "after_mar10_treatSES3", "after_mar10_treatSES4", "after_mar10_treatSES5")))
Wealth<-as.matrix(subset(data, select=c("SES1", "SES2", "SES3", "SES4", "SES5")))
FedeligW<-as.matrix(subset(data, select=c("fedeligSES1", "fedeligSES2", "fedeligSES3", "fedeligSES4", "fedeligSES5")))


####################################################################################
# CLUSTERED ROBUST
####################################################################################
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  Robust<-as.matrix(coeftest(fm, vcovCL))
  return(Robust) }


cl2 <- function(dat,fm, cluster){
  attach(dat, warn.conflicts = F)
  require(sandwich)
  require(lmtest)
  not <- attr(fm$model,"na.action")
  cluster <- cluster[-not]
  with( dat[-not,] ,{
    M <- length(unique(cluster))
    N <- length(cluster)
    K <- fm$rank
    dfc <- (M/(M-1))*((N-1)/(N-K))
    uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
    vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
    coeftest(fm, vcovCL)
  }
  )
}


####################################################################################
#Month level charts
####################################################################################

ExpVarsNoTime  <- as.matrix(subset(data, select=c(age17,age18,age19,age20,age21,age22,age23,age24,age25,age27,age28,age29,female,hispanic,white,asian,other,student,mar,ue,fpl_ratio,fpl_ratio_2,ue_treat)))
EffectsN<-as.matrix(subset(data, select=c(mar_sept10, after_oct10,elig_mar10,elig_oct10)))
MonthYear <- as.matrix(paste(data$year,data$monthchar,sep=","))

# Run original reg, no subgroups
Results<-cl(data,lm(Outcomes[,1]~ExpVarsNoTime+factor(data$monthdate)+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table11a <- rbind(c(0,0,0,0),as.matrix(Results)[25:63,1:4])
table11a[,1]<-table11a[,1]+Results[1,1]
Results<-cl(data,lm(Outcomes[,2]~ExpVarsNoTime+factor(data$monthdate)+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table11b <- rbind(c(0,0,0,0),as.matrix(Results)[25:63,1:4])
table11b[,1]<-table11b[,1]+Results[1,1]
Results<-cl(data,lm(Outcomes[,3]~ExpVarsNoTime+factor(data$monthdate)+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table11c <- rbind(c(0,0,0,0),as.matrix(Results)[25:63,1:4])
table11c[,1]<-table11c[,1]+Results[1,1]
Results<-cl(data,lm(Outcomes[,4]~ExpVarsNoTime+factor(data$monthdate)+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table11d <- rbind(c(0,0,0,0),as.matrix(Results)[25:63,1:4])
table11d[,1]<-table11d[,1]+Results[1,1]
Results<-cl(data,lm(Outcomes[,5]~ExpVarsNoTime+factor(data$monthdate)+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table11e <- rbind(c(0,0,0,0),as.matrix(Results)[25:63,1:4])
table11e[,1]<-table11e[,1]+Results[1,1]
dates <- as.data.frame(c(as.Date("2008-8-1"),as.Date(substr(rownames(table11a[2:nrow(table11a),]), 23, 32))))

colnames(dates) <- "dates"
table11a<-cbind(dates,as.data.frame(table11a))
table11b<-cbind(dates,as.data.frame(table11b))
table11c<-cbind(dates,as.data.frame(table11c))
table11d<-cbind(dates,as.data.frame(table11d))
table11e<-cbind(dates,as.data.frame(table11e))
colnames(table11a)<-colnames(table11b)<-colnames(table11c)<-colnames(table11d)<-colnames(table11e)<-c("dates","Estimate","SE","tvalue","pvalue")

p11a <- ggplot(table11a, aes(dates)) + ggtitle("Any source") +
  geom_line(aes(y=Estimate)) + 
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") + 
  scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months")) +
  theme(axis.ticks = element_blank(), axis.text.x = element_blank(),axis.title.x=element_blank()) + ylab("intercept + month dummy")
p11b <- ggplot(table11b, aes(dates)) + ggtitle("Employer dependent coverage (through parents)") +
  geom_line(aes(y=Estimate)) + 
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") +
  scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months")) + 
  theme(axis.ticks = element_blank(), axis.text.x = element_blank(),axis.title.x=element_blank()) + ylab("intercept + month dummy")
p11c <- ggplot(table11c, aes(dates)) + ggtitle("Individually purchased insurance in own name") +
  geom_line(aes(y=Estimate)) + 
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") +
  scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months")) + 
  theme(axis.ticks = element_blank(), axis.text.x = element_blank(),axis.title.x=element_blank()) + ylab("intercept + month dummy")
p11d <- ggplot(table11d, aes(dates)) + ggtitle("Employer own coverage") +
  geom_line(aes(y=Estimate)) + 
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") +
  scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months")) + 
  theme(axis.ticks = element_blank(), axis.text.x = element_blank(),axis.title.x=element_blank()) + ylab("intercept + month dummy")
p11e <- ggplot(table11e, aes(dates)) + ggtitle("Government provided") +
  geom_line(aes(y=Estimate)) + 
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") +
  scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months")) + 
  xlab ("date") + ylab("intercept + month dummy")

jpeg("ModTables/time trends.jpg",width = 800, height = 1000, quality = 100)
grid.arrange(p11a, p11b,p11c,p11d,p11e, nrow=5)
dev.off()

# Trend + level shifts
Results<-cl(data,lm(Outcomes[,1]~ExpVarsNoTime+data$trend+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table12a <- as.matrix(Results)[c(1:1,25),1:4]
Results<-cl(data,lm(Outcomes[,1]~ExpVarsNoTime+data$trend+data$trend2+data$fedelig+factor(data$fipstate), w=data$p_weight), data$fipstate)
table13a <- as.matrix(Results)[c(1:1,25:26),1:4]

# Predicted fit
pred <- as.data.frame(matrix(0,nrow=max(data$trend)+1,ncol=4))
colnames(pred) <- c("date","pred_lin","pred_quad","pred_spline")
pred$date<-as.Date("2000-1-1")
for(i in 0:max(data$trend)){
  imonth1<-8+i
  iyear<-2008+trunc((imonth1-1)/12)
  imonth<-8+i - (iyear-2008)*12
  idate <- as.Date(paste(iyear,imonth,"1",sep="-"))
  t_pred_lin <- table12a[1,1]
  t_pred_lin<-t_pred_lin+i*table12a[2,1]
  t_pred_quad <- table13a[1,1]
  t_pred_quad<-t_pred_quad+i*table13a[2,1]+i*i*table13a[3,1]
  
  pred$pred_lin[i+1] <- t_pred_lin
  pred$pred_quad[i+1] <- t_pred_quad
  pred$date[i+1] <- idate
}
trend_assumps <- as.data.frame(cbind(table11a,pred$pred_lin,pred$pred_quad))
colnames(trend_assumps) <- c(colnames(table11a),"pred_lin","pred_quad")

p12a <- ggplot(trend_assumps, aes(dates)) + ggtitle("Any source, Trend Assumptions") +
  geom_line(aes(y=Estimate)) +
  geom_ribbon(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), alpha=0.2,fill="blue") +
  geom_line(aes(y=pred_lin)) + geom_line(aes(y=pred_quad))
scale_x_date(breaks = date_breaks("years"), minor_breaks = date_breaks("months"))


jpeg("ModTables/trend assumptions.jpg", quality = 100)
p12a
dev.off()

####################################################################################
# Time trends, original variables of interest
####################################################################################

TTSens<-matrix(NA,6,8)
colnames(TTSens)<-c("Any Source, linear trend", "Employer dependent coverage (through parents), linear trend","Any Source, quadratic trend", "Employer dependent coverage (through parents), quadratic trend","Any Source, state linear trend", "Employer dependent coverage (through parents), state linear trend","Any Source, state quadratic trend", "Employer dependent coverage (through parents), state quadratic trend")
rownames(TTSens)<-c("ACA enactment effect (March-Sept 2011)", "ACA implementation effect (October 2010-Novermber 2011)", "Treatment, before ACA enactment", "Control, before", "Treatment, after ACA implementation", "Control, after")

# Linear trend
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~Effects+ExpVarsNoTime+data$trend+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[5,4]<=0.01) stars1<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars1<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[6,4]<=0.01) stars2<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars2<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars2<-"*" else stars2<-""
  
  TTSens[1,i]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars1, sep="") 
  TTSens[2,i]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars2, sep="")                   
  
}
rm(i, Results, stars1, stars2)
TTSens[3:6,1]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$anyhi,x$p_weight)))[1:4,"Res"], digits=3)
TTSens[3:6,2]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$emphi_dep,x$p_weight)))[1:4,"Res"], digits=3)

# Quadratic trend
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~Effects+ExpVarsNoTime+data$trend+data$trend2+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[5,4]<=0.01) stars1<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars1<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[6,4]<=0.01) stars2<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars2<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars2<-"*" else stars2<-""
  
  TTSens[1,i+2]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars1, sep="") 
  TTSens[2,i+2]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars2, sep="")                   
  
}
rm(i, Results, stars1, stars2)
TTSens[3:6,3]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$anyhi,x$p_weight)))[1:4,"Res"], digits=3)
TTSens[3:6,4]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$emphi_dep,x$p_weight)))[1:4,"Res"], digits=3)

# State linear trend
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~Effects+ExpVarsNoTime+StateTrends+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[5,4]<=0.01) stars1<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars1<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[6,4]<=0.01) stars2<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars2<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars2<-"*" else stars2<-""
  
  TTSens[1,i+4]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars1, sep="") 
  TTSens[2,i+4]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars2, sep="")                   
  
}
rm(i, Results, stars1, stars2)
TTSens[3:6,5]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$anyhi,x$p_weight)))[1:4,"Res"], digits=3)
TTSens[3:6,6]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$emphi_dep,x$p_weight)))[1:4,"Res"], digits=3)

# State linear trend
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~Effects+ExpVarsNoTime+StateTrends+StateTrend2s+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[5,4]<=0.01) stars1<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars1<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[6,4]<=0.01) stars2<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars2<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars2<-"*" else stars2<-""
  
  TTSens[1,i+6]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars1, sep="") 
  TTSens[2,i+6]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars2, sep="")                   
  
}
rm(i, Results, stars1, stars2)
TTSens[3:6,7]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$anyhi,x$p_weight)))[1:4,"Res"], digits=3)
TTSens[3:6,8]<-round(ddply(data, .(sum_group), function(x) data.frame(Res=weighted.mean(x$emphi_dep,x$p_weight)))[1:4,"Res"], digits=3)

TTSens
write.csv(file="ModTables/TTSens.csv", x=TTSens)
print.xtable(TTSens, type="latex",file = "ModTables/TTSens.tex")


####################################################################################
# Time trends, alternative specifications
####################################################################################

TimeSpecs<-matrix("",7,12)
#colnames(TimeSpecs)<-c("Any Source, levels", "Employer dependent coverage (through parents), levels","Any Source, levels and state trends", "Employer dependent coverage (through parents), levels and state trends","Any Source, spline and state trends", "Employer dependent coverage (through parents), spline and state trends","Any Source, levels, spline, state trends", "Employer dependent coverage (through parents), levels, spline, state trends","Any Source, state quadratic trend", "Employer dependent coverage (through parents), state quadratic trend")
rownames(TimeSpecs)<-c("Eligible","After March 10", "After March 10 * Eligible", "Trend after March 10", "Trend after March 10 * Eligible","State linear trends","State linear trends^2")

# No trend
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+EffectsSplineLevels+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  
  TimeSpecs[1,i]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[2,i]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[3,i]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")                   
  TimeSpecs[6,i]<-"No"
  TimeSpecs[7,i]<-"No"
}
rm(i, Results, stars1, stars2,stars3)

# Levels and state trends
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+EffectsSplineLevels+StateTrends+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  
  TimeSpecs[1,i+2]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[2,i+2]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[3,i+2]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")                   
  TimeSpecs[6,i+2]<-"Yes"
  TimeSpecs[7,i+2]<-"No"
}
rm(i, Results, stars1, stars2,stars3)

# Linear trend break, state trends
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+data$trendafter_mar10+data$trendafter_mar10_treat+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  
  TimeSpecs[1,i+4]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[4,i+4]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[5,i+4]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")                   
  TimeSpecs[6,i+4]<-"No"
  TimeSpecs[7,i+4]<-"No"
}
rm(i, Results, stars1, stars2,stars3)


# Linear trend break, state trends
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+data$trendafter_mar10+data$trendafter_mar10_treat+StateTrends+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-as.character(NA)
  stars2<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  
  TimeSpecs[1,i+6]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[4,i+6]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[5,i+6]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")                   
  TimeSpecs[6,i+6]<-"Yes"
  TimeSpecs[7,i+6]<-"No"
}
rm(i, Results, stars1, stars2,stars3)

# Linear trend break, state trends
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+EffectsSplineLevels+data$trendafter_mar10+data$trendafter_mar10_treat+StateTrends+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-stars4<-stars5<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  
  TimeSpecs[1,i+8]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[2,i+8]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[3,i+8]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")  
  TimeSpecs[4,i+8]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars4, sep="")                   
  TimeSpecs[5,i+8]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars5, sep="")                   
  TimeSpecs[6,i+8]<-"Yes"
  TimeSpecs[7,i+8]<-"No"
}
rm(i, Results, stars1, stars2,stars3,stars4,stars5)

# Linear trend break, state trends^2
for(i in 1:2){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+EffectsSplineLevels+data$trendafter_mar10+data$trendafter_mar10_treat+StateTrends+StateTrend2s+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  stars1<-stars2<-stars3<-stars4<-stars5<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  
  TimeSpecs[1,i+10]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  TimeSpecs[2,i+10]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  TimeSpecs[3,i+10]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")  
  TimeSpecs[4,i+10]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars4, sep="")                   
  TimeSpecs[5,i+10]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars5, sep="")                   
  TimeSpecs[6,i+10]<-"Yes"
  TimeSpecs[7,i+10]<-"Yes"
}
rm(i, Results, stars1, stars2,stars3,stars4,stars5)

# Format for output

TimeSpecsOut <- matrix("",nrow(TimeSpecs),ncol(TimeSpecs))
for (i in 1:ncol(TimeSpecs)/2){
  TimeSpecsOut[,i] <- TimeSpecs[,i*2-1]
  TimeSpecsOut[,i+ncol(TimeSpecs)/2] <- TimeSpecs[,i*2]
}

TimeSpecsOut <- cbind(rownames(TimeSpecs),TimeSpecsOut)
TimeSpecsOutNames <- t(as.matrix(c("","Any Source, levels","","","","","", "Employer dependent coverage (through parents), levels","","","","","")))
TimeSpecsOut <- rbind(TimeSpecsOutNames,TimeSpecsOut)
rownames(TimeSpecsOut) <- TimeSpecsOut[,1]

write.csv(file="ModTables/TimeSpecs.csv", x=TimeSpecsOut)
print.xtable(TimeSpecsOut, type="latex",file = "ModTables/TimeSpecs.tex")




####################################################################################
# Table 1 Descriptive Statistics
####################################################################################

Table1Var<-c("anyhi","emphi_dep","emphi","indiv", "govhi", "employed", "unemployed", "age", "white", "black", "hispanic", "mar", "ft_student","hsdo","hsg", "somcol", "colgrd", "bad_hlth", "famyincy")
Names<-c("Covered by any health insurance (HI)", "Covered as parent's dependent", "Covered by own employer HI", "Covered by individually purchased HI in own name", "Covered by government HI", "Employed", "Unemployed","Age", "White", "African-American", "Hispanic", "Married", "Student", "Less than high school", "High school graduate", "Some college", "College graduate", "Self-reported health is less than excellent", "Average family yearly income", "Number of person-month observations", "Corresponding number of unique persons")

Table1<-matrix(0,length(Names), 4)
rownames(Table1)<-Names
colnames(Table1)<-c("All", "Age 16-18", "Age 19-25", "Age 27-29")

for(i in 1: length(Table1Var)){
  
  Table1[i,1]<-round(weighted.mean(data[,Table1Var[i]],data$p_weight, na.rm=T),digits=3)
  Table1[i,2]<-round(weighted.mean(data[data$age_group==1,Table1Var[i]],data[data$age_group==1,]$p_weight, na.rm=T),digits=3)
  Table1[i,3]<-round(weighted.mean(data[data$age_group==2,Table1Var[i]],data[data$age_group==2,]$p_weight, na.rm=T),digits=3)
  Table1[i,4]<-round(weighted.mean(data[data$age_group==3,Table1Var[i]],data[data$age_group==3,]$p_weight, na.rm=T),digits=3)
}

Table1[20,1]<-round(nrow(data), digits=0)
Table1[20,2]<-round(sum(as.integer(data$age_group==1)), digits=0)
Table1[20,3]<-round(sum(as.integer(data$age_group==2)), digits=0)
Table1[20,4]<-round(sum(as.integer(data$age_group==3)), digits=0)
Table1[21,1]<-round(length(unique(data$groupid)), digits=0)
Table1[21,2]<-round(length(unique(data[data$age_group==1,]$groupid)), digits=0)
Table1[21,3]<-round(length(unique(data[data$age_group==2,]$groupid)), digits=0)
Table1[21,4]<-round(length(unique(data[data$age_group==3,]$groupid)), digits=0)

rm(i, Names, Table1Var)
Table1

write.csv(file="ModTables/Descriptive.csv", x=Table1)
Table1<- xtable(Table1)
print.xtable(Table1, type="latex",file = "ModTables/Descriptive.tex")

####################################################################################
# Table 2 Descriptive Statistics
####################################################################################
Var<-c("anyhi","emphi_dep","indiv","emphi", "govhi")
Table2<-matrix("",16,5)
colnames(Table2)<-c("Any Source", "Employer dependent coverage (through parents)", "Individually purchased insurance in own name", "Employer own coverage", "Government provided")
rownames(Table2)<-c("Health: Excellent", "Very Good", "Good", "Fair", "Poor", "Wealth: 1st Quintile [Lowest]", "2nd Quintile", "3rd Quintile", "4th Quintile", "Wealth: 5th Quintile [Highest]", "Employed", "Full Time", "Part-Time", "Variable Hours", "Change Job Status","Change Employer")

for(i in 1: length(Var)){
  
  Table2[1,i]<-round(weighted.mean(data[data$ExH==1,Var[i]],data[data$ExH==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[2,i]<-round(weighted.mean(data[data$VGH==1,Var[i]],data[data$VGH==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[3,i]<-round(weighted.mean(data[data$GH==1 ,Var[i]],data[data$GH==1 , "p_weight"], na.rm=T),digits=3)
  Table2[4,i]<-round(weighted.mean(data[data$FH==1 ,Var[i]],data[data$FH==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[5,i]<-round(weighted.mean(data[data$PH==1 ,Var[i]],data[data$PH==1 ,"p_weight"], na.rm=T),digits=3)
  
  Table2[6,i]<-round(weighted.mean(data[data$SES1==1 ,Var[i]],data[data$SES1==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[7,i]<-round(weighted.mean(data[data$SES2==1 ,Var[i]],data[data$SES2==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[8,i]<-round(weighted.mean(data[data$SES3==1 ,Var[i]],data[data$SES3==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[9,i]<-round(weighted.mean(data[data$SES4==1 ,Var[i]],data[data$SES4==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[10,i]<-round(weighted.mean(data[data$SES5==1 ,Var[i]],data[data$SES5==1,"p_weight"], na.rm=T),digits=3)
  
  Table2[11,i]<-round(weighted.mean(data[data$employed==1 ,Var[i]],data[data$employed==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[12,i]<-round(weighted.mean(data[data$ft==1 ,Var[i]],data[data$ft==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[13,i]<-round(weighted.mean(data[data$pt==1 ,Var[i]],data[data$pt==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[14,i]<-round(weighted.mean(data[data$hour_vary==1 ,Var[i]],data[data$hour_vary==1,"p_weight"], na.rm=T),digits=3)
  Table2[15,i]<-round(weighted.mean(data[data$job_ease==1 ,Var[i]],data[data$job_ease==1 ,"p_weight"], na.rm=T),digits=3)
  Table2[16,i]<-round(weighted.mean(data[data$emp2diffemp==1 ,Var[i]],data[data$emp2diffemp==1 ,"p_weight"], na.rm=T),digits=3)
    
}


rm(i, Var)

write.csv(file="ModTables/DescriptiveSelection.csv", x=Table2)
Table2<- xtable(Table2)
print.xtable(Table2, type="latex",file = "ModTables/DescriptiveSelection.tex")


####################################################################################
# Table 3 Results with Modified Model
####################################################################################

Table3<-matrix("",7,5)
colnames(Table3)<-c("Any Source", "Employer dependent coverage (through parents)", "Individually purchased insurance in own name", "Employer own coverage", "Government provided")
rownames(Table3)<-c("Eligible","After March 10", "After Mar'10*Eligible", "Trend after Mar'10", "Trend after Mar'10*Eligible","State linear trends","State linear trends^2")

for(i in 1:5){  
  Results<-cl(data,lm(Outcomes[,i]~data$fedelig+EffectsSplineLevels+data$trendafter_mar10+data$trendafter_mar10_treat+StateTrends+StateTrend2s+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate)
  
  stars1<-stars2<-stars3<-stars4<-stars5<-as.character(NA)
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  
  Table3[1,i]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="") 
  Table3[2,i]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="")                   
  Table3[3,i]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")  
  Table3[4,i]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars4, sep="")                   
  Table3[5,i]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars5, sep="")                   
  Table3[6,i]<-"Yes"
  Table3[7,i]<-"Yes"
}

rm(i, Results, stars1, stars2,stars3,stars4,stars5)
Table3

write.csv(file="ModTables/General Results Table.csv", x=Table3)
Table3<- xtable(Table3)
print.xtable(Table3, type="latex",file = "ModTables/General Results Table.tex")

####################################################################################
# Table 4 Results for Health Selection
####################################################################################
Table4<-matrix(NA,12,5)
colnames(Table4)<-c("Any Source", "Employer dependent coverage (through parents)", "Individually purchased insurance in own name", "Employer own coverage", "Government provided")
rownames(Table4)<-c("Trend after March 10: Excellent Health", "Very Good Health", "Good Health", "Fair Health", "Poor Health", "Trend after March 10*Eligible: Excellent Health", "Very Good Health", "Good Health", "Fair Health", "Poor Health", "State linear trends","State linear trends^2")

for(i in 1:ncol(Outcomes)){
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  stars3<-as.character(NA)
  stars4<-as.character(NA)
  stars5<-as.character(NA)
  stars6<-as.character(NA)
  stars7<-as.character(NA)
  stars8<-as.character(NA)
  stars9<-as.character(NA)
  stars10<-as.character(NA)
  
  Results<-cl2(data,lm(Outcomes[,i]~EffectsH+FedeligH+EffectsSplineLevelsH+Health+StateTrends+StateTrend2s+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate) 
  
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  if (Results[7,4]<=0.01) stars6<-"***" else if (Results[7,4]<=0.05 & Results[7,4]>0.01) stars6<-"**" else if (Results[7,4]<=0.10 & Results[7,4]>0.05) stars6<-"*" else stars6<-""
  if (Results[8,4]<=0.01) stars7<-"***" else if (Results[8,4]<=0.05 & Results[8,4]>0.01) stars7<-"**" else if (Results[8,4]<=0.10 & Results[8,4]>0.05) stars7<-"*" else stars7<-""
  if (Results[9,4]<=0.01) stars8<-"***" else if (Results[9,4]<=0.05 & Results[9,4]>0.01) stars8<-"**" else if (Results[9,4]<=0.10 & Results[9,4]>0.05) stars8<-"*" else stars8<-""
  if (Results[10,4]<=0.01) stars9<-"***" else if (Results[10,4]<=0.05 & Results[10,4]>0.01) stars9<-"**" else if (Results[10,4]<=0.10 & Results[10,4]>0.05) stars9<-"*" else stars9<-""
  if (Results[11,4]<=0.01) stars10<-"***" else if (Results[11,4]<=0.05 & Results[11,4]>0.01) stars10<-"**" else if (Results[11,4]<=0.10 & Results[11,4]>0.05) stars10<-"*" else stars10<-""
  
  
  Table4[1,i]<-paste(format(round(Results[2,1], digits=4),scientific=FALSE), " ","(",round(Results[2,2], digits=4),")", stars1, sep="")
  Table4[2,i]<-paste(format(round(Results[3,1], digits=4),scientific=FALSE), " ","(",round(Results[3,2], digits=4),")", stars2, sep="") 
  Table4[3,i]<-paste(format(round(Results[4,1], digits=4),scientific=FALSE), " ","(",round(Results[4,2], digits=4),")", stars3, sep="")
  Table4[4,i]<-paste(format(round(Results[5,1], digits=4),scientific=FALSE), " ","(",round(Results[5,2], digits=4),")", stars4, sep="")
  Table4[5,i]<-paste(format(round(Results[6,1], digits=4),scientific=FALSE), " ","(",round(Results[6,2], digits=4),")", stars5, sep="")
  Table4[6,i]<-paste(format(round(Results[7,1], digits=4),scientific=FALSE), " ","(",round(Results[7,2], digits=4),")", stars6, sep="") 
  Table4[7,i]<-paste(format(round(Results[8,1], digits=4),scientific=FALSE), " ","(",round(Results[8,2], digits=4),")", stars7, sep="")
  Table4[8,i]<-paste(format(round(Results[9,1], digits=4),scientific=FALSE), " ","(",round(Results[9,2], digits=4),")", stars8, sep="") 
  Table4[9,i]<-paste(format(round(Results[10,1], digits=4),scientific=FALSE), " ","(",round(Results[10,2], digits=4),")", stars9, sep="")
  Table4[10,i]<-paste(format(round(Results[11,1], digits=4),scientific=FALSE), " ","(",round(Results[11,2], digits=4),")", stars10, sep="")
  Table4[11,i]<-"Yes"
  Table4[12,i]<-"Yes"
  
}

Table4

write.csv(file="ModTables/Health Selection Results.csv", x=Table4)
Table4<- xtable(Table4)
print.xtable(Table4, type="latex",file = "ModTables/Health Selection Results.tex")

rm(i, Results, stars1, stars2, stars3, stars4, stars5, stars6, stars7, stars8, stars9, stars10, EffectsH, EffectsSplineLevelsH, FedeligH, Health)

####################################################################################
# Table 5 Results for Wealth Selection
####################################################################################
Table5<-matrix(NA,12,5)
colnames(Table5)<-c("Any Source", "Employer dependent coverage (through parents)", "Individually purchased insurance in own name", "Employer own coverage", "Government provided")
rownames(Table5)<-c("Trend after March 10: 1st Quintile [Lowest]", "Trend after March 10: 2nd Quintile", "Trend after March 10: 3rd Quintile", "Trend after March 10: 4th Quintile", "Trend after March 10: 5th Quintile [Highest]", "Trend after March 10 * Eligible: 1st Quintile [Lowest]", "Trend after March 10 * Eligible: 2nd Quintile", "Trend after March 10 * Eligible: 3rd Quintile", "Trend after March 10 * Eligible: 4rd Quintile", "Trend after March 10 * Eligible: 5th Quintile [Highest]", "State linear trends","State linear trends^2")

for(i in 1:ncol(Outcomes)){
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  stars3<-as.character(NA)
  stars4<-as.character(NA)
  stars5<-as.character(NA)
  stars6<-as.character(NA)
  stars7<-as.character(NA)
  stars8<-as.character(NA)
  stars9<-as.character(NA)
  stars10<-as.character(NA)
  
  Results<-cl(data,lm(Outcomes[,i]~EffectsW+FedeligW+EffectsSplineLevelsW+Wealth+StateTrends+StateTrend2s+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate) 
  
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  if (Results[7,4]<=0.01) stars6<-"***" else if (Results[7,4]<=0.05 & Results[7,4]>0.01) stars6<-"**" else if (Results[7,4]<=0.10 & Results[7,4]>0.05) stars6<-"*" else stars6<-""
  if (Results[8,4]<=0.01) stars7<-"***" else if (Results[8,4]<=0.05 & Results[8,4]>0.01) stars7<-"**" else if (Results[8,4]<=0.10 & Results[8,4]>0.05) stars7<-"*" else stars7<-""
  if (Results[9,4]<=0.01) stars8<-"***" else if (Results[9,4]<=0.05 & Results[9,4]>0.01) stars8<-"**" else if (Results[9,4]<=0.10 & Results[9,4]>0.05) stars8<-"*" else stars8<-""
  if (Results[10,4]<=0.01) stars9<-"***" else if (Results[10,4]<=0.05 & Results[10,4]>0.01) stars9<-"**" else if (Results[10,4]<=0.10 & Results[10,4]>0.05) stars9<-"*" else stars9<-""
  if (Results[11,4]<=0.01) stars10<-"***" else if (Results[11,4]<=0.05 & Results[11,4]>0.01) stars10<-"**" else if (Results[11,4]<=0.10 & Results[11,4]>0.05) stars10<-"*" else stars10<-""
  
  
  Table5[1,i]<-paste(round(Results[2,1], digits=4), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="")
  Table5[2,i]<-paste(round(Results[3,1], digits=4), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="") 
  Table5[3,i]<-paste(round(Results[4,1], digits=4), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")
  Table5[4,i]<-paste(round(Results[5,1], digits=4), " ","(",(round(Results[5,2], digits=4)),")", stars4, sep="")
  Table5[5,i]<-paste(round(Results[6,1], digits=4), " ","(",(round(Results[6,2], digits=4)),")", stars5, sep="")
  Table5[6,i]<-paste(round(Results[7,1], digits=4), " ","(",(round(Results[7,2], digits=4)),")", stars6, sep="") 
  Table5[7,i]<-paste(round(Results[8,1], digits=4), " ","(",(round(Results[8,2], digits=4)),")", stars7, sep="")
  Table5[8,i]<-paste(round(Results[9,1], digits=4), " ","(",(round(Results[9,2], digits=4)),")", stars8, sep="") 
  Table5[9,i]<-paste(round(Results[10,1], digits=4), " ","(",(round(Results[10,2], digits=4)),")", stars9, sep="")
  Table5[10,i]<-paste(round(Results[11,1], digits=4), " ","(",(round(Results[11,2], digits=4)),")", stars10, sep="")
  Table5[11,i]<-"Yes"
  Table5[12,i]<-"Yes"
  
}

Table5

write.csv(file="ModTables/Wealth Selection Results.csv", x=Table5)
Table5<- xtable(Table5)
print.xtable(Table5, type="latex",file = "Wealth Selection Results.tex")

rm(i, Results, stars1, stars2, stars3, stars4, stars5, stars6, stars7, stars8, stars9, stars10, EffectsW, EffectsSplineLevelsW, FedeligW, Wealth)


####################################################################################
#Selection Issues with Employment
####################################################################################



Table5<-matrix(NA,14,5)
colnames(Table5)<-c("Any Source", "Employer dependent coverage (through parents)", "Individually purchased insurance in own name", "Employer own coverage", "Government provided")
rownames(Table5)<-c("Trend after March 10: Employed", "Full Time", "Part Time", "Variable Hours", "Change Job Status", "Change Employer", "Trend after March 10 * Eligible: Employed", "*Full Time","*Part Time", "*Variable Hours", "*Change Job Status", "*Change Employer", "State linear trends","State linear trends^2")

for(i in 1:ncol(Outcomes)){
  stars1<-as.character(NA)
  stars2<-as.character(NA)
  stars3<-as.character(NA)
  stars4<-as.character(NA)
  stars5<-as.character(NA)
  stars6<-as.character(NA)
  stars7<-as.character(NA)
  stars8<-as.character(NA)
  stars9<-as.character(NA)
  stars10<-as.character(NA)
  stars11<-as.character(NA)
  stars12<-as.character(NA)
  Results<-cl2(data,lm(Outcomes[,i]~EffectsE+FedeligE+EffectsSplineLevelsE+employed+StateTrends+StateTrend2s+ExpVarsNoTime+factor(data$fipstate), w=data$p_weight), data$fipstate) 
  
  if (Results[2,4]<=0.01) stars1<-"***" else if (Results[2,4]<=0.05 & Results[2,4]>0.01) stars1<-"**" else if (Results[2,4]<=0.10 & Results[2,4]>0.05) stars1<-"*" else stars1<-""
  if (Results[3,4]<=0.01) stars2<-"***" else if (Results[3,4]<=0.05 & Results[3,4]>0.01) stars2<-"**" else if (Results[3,4]<=0.10 & Results[3,4]>0.05) stars2<-"*" else stars2<-""
  if (Results[4,4]<=0.01) stars3<-"***" else if (Results[4,4]<=0.05 & Results[4,4]>0.01) stars3<-"**" else if (Results[4,4]<=0.10 & Results[4,4]>0.05) stars3<-"*" else stars3<-""
  if (Results[5,4]<=0.01) stars4<-"***" else if (Results[5,4]<=0.05 & Results[5,4]>0.01) stars4<-"**" else if (Results[5,4]<=0.10 & Results[5,4]>0.05) stars4<-"*" else stars4<-""
  if (Results[6,4]<=0.01) stars5<-"***" else if (Results[6,4]<=0.05 & Results[6,4]>0.01) stars5<-"**" else if (Results[6,4]<=0.10 & Results[6,4]>0.05) stars5<-"*" else stars5<-""
  if (Results[7,4]<=0.01) stars6<-"***" else if (Results[7,4]<=0.05 & Results[7,4]>0.01) stars6<-"**" else if (Results[7,4]<=0.10 & Results[7,4]>0.05) stars6<-"*" else stars6<-""
  if (Results[8,4]<=0.01) stars7<-"***" else if (Results[8,4]<=0.05 & Results[8,4]>0.01) stars7<-"**" else if (Results[8,4]<=0.10 & Results[8,4]>0.05) stars7<-"*" else stars7<-""
  if (Results[9,4]<=0.01) stars8<-"***" else if (Results[9,4]<=0.05 & Results[9,4]>0.01) stars8<-"**" else if (Results[9,4]<=0.10 & Results[9,4]>0.05) stars8<-"*" else stars8<-""
  if (Results[10,4]<=0.01) stars9<-"***" else if (Results[10,4]<=0.05 & Results[10,4]>0.01) stars9<-"**" else if (Results[10,4]<=0.10 & Results[10,4]>0.05) stars9<-"*" else stars9<-""
  if (Results[11,4]<=0.01) stars10<-"***" else if (Results[11,4]<=0.05 & Results[11,4]>0.01) stars10<-"**" else if (Results[11,4]<=0.10 & Results[11,4]>0.05) stars10<-"*" else stars10<-""
  if (Results[12,4]<=0.01) stars11<-"***" else if (Results[12,4]<=0.05 & Results[12,4]>0.01) stars11<-"**" else if (Results[12,4]<=0.10 & Results[12,4]>0.05) stars11<-"*" else stars11<-""
  if (Results[13,4]<=0.01) stars12<-"***" else if (Results[13,4]<=0.05 & Results[13,4]>0.01) stars12<-"**" else if (Results[13,4]<=0.10 & Results[13,4]>0.05) stars2<-"*" else stars12<-""
  
  
  Table5[1,i]<-paste(format(round(Results[2,1], digits=4),scientific=FALSE), " ","(",(round(Results[2,2], digits=4)),")", stars1, sep="")
  Table5[2,i]<-paste(format(round(Results[3,1], digits=4),scientific=FALSE), " ","(",(round(Results[3,2], digits=4)),")", stars2, sep="") 
  Table5[3,i]<-paste(format(round(Results[4,1], digits=4),scientific=FALSE), " ","(",(round(Results[4,2], digits=4)),")", stars3, sep="")
  Table5[4,i]<-paste(format(round(Results[5,1], digits=4),scientific=FALSE), " ","(",(round(Results[5,2], digits=4)),")", stars4, sep="")
  Table5[5,i]<-paste(format(round(Results[6,1], digits=4),scientific=FALSE), " ","(",(round(Results[6,2], digits=4)),")", stars5, sep="")
  Table5[6,i]<-paste(format(round(Results[7,1], digits=4),scientific=FALSE), " ","(",(round(Results[7,2], digits=4)),")", stars6, sep="") 
  Table5[7,i]<-paste(format(round(Results[8,1], digits=4),scientific=FALSE), " ","(",(round(Results[8,2], digits=4)),")", stars7, sep="")
  Table5[8,i]<-paste(format(round(Results[9,1], digits=4),scientific=FALSE), " ","(",(round(Results[9,2], digits=4)),")", stars8, sep="") 
  Table5[9,i]<-paste(format(round(Results[10,1], digits=4),scientific=FALSE), " ","(",(round(Results[10,2], digits=4)),")", stars9, sep="")
  Table5[10,i]<-paste(format(round(Results[11,1], digits=4),scientific=FALSE), " ","(",(round(Results[11,2], digits=4)),")", stars10, sep="")
  Table5[11,i]<-paste(format(round(Results[12,1], digits=4),scientific=FALSE), " ","(",(round(Results[12,2], digits=4)),")", stars11, sep="")
  Table5[12,i]<-paste(format(round(Results[13,1], digits=4),scientific=FALSE), " ","(",(round(Results[13,2], digits=4)),")", stars12, sep="")
  Table5[13,i]<-"Yes"
  Table5[14,i]<-"Yes"
  
}

Table5

write.csv(file="Employement Selection Results Table 5.csv", x=Table5)
Table5<- xtable(Table5)
print.xtable(Table5, type="latex",file = "Employment Selection Results Table 5.tex")


